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This contribution reports on numerical simulations of 2D granular flows on erodible beds. The 
broad aim is to investigate whether simple flows of model granular matter exhibits spontaneous 
oscillatory motion in generic flow conditions, and in this case, whether the frictional properties of 
the contacts between grains may affect the existence or the characteristics of this oscillatory motion. 
The analysis of different series of simulations show that the flow develops an oscillatory motion with 
a well-defined frequency which increases like the inverse of the velocity's square root. We show that 
the oscillation is essentially a surface phenomena. The amplitude of the oscillation is higher for 
lower volume fractions, and can thus be related to the flow velocity and grains friction properties. 
The study of the influence of the periodic geometry of the simulation cell shows no significant effect. 
These results are discussed in relation to sonic sands. 

PACS numbers: 45.70.-n, 05.65.+b 



I. INTRODUCTION 



Granular flows have long been the subject of sustained 
interest due to their surprising properties: the existence 
of an internal yield stress allowing them both solid-like 
and fluid-like behaviors and more generally their elu- 
sive rheology [TH2] , their ability to segregate according to 
grain size or mass [T0ifT4] , their similarities with glassy 
systems [T5HT8] . their non-local properties [19-22 j, to cite 
only few examples among many. While most of the fea- 
tures listed above have been the subject of careful lab- 
oratory experiments, some peculiarities of granular sys- 
tems express themselves best in nature: this is the case 
of musical (or sonic) sands, and most famously, of boom- 
ing dunes [23-29 . Booming dunes, when a surface flow 
is triggered at their flank, can emit a surprisingly loud 
sound which has proven a long-lasting challenge to ex- 
perimentalist and theoretician alike. In spite of numerous 
advances, including experimental investigation, modeling 
and field measurements, the mechanisms at play in the 
booming phenomena remain controversial (see the recent 
review by [30] and reference therein). However, converg- 
ing observations have emerged. It was observed that the 
granular surface flow alone can produce the sound, with- 
out the dune being a necessary ingredient. Then, both 
laboratory experiments and field measurement show that 
a minimum flow velocity is needed for sound to be emit- 
ted. Moreover, sonic sand grains need specific contact 
properties, ensured in nature by a mineral coating repro- 
duced in the lab by successive bath in complex salty so- 
lutions [28]. Based on these observations, two questions 
arises: can a simple granular flow spontaneously exhibit 
oscillatory motion with a well-defined frequency? If yes, 
how do frictional contact properties affect the oscillatory 
dynamics? These two questions are the subject of the 
present contribution. While intermittent motion close to 
jamming transition was already observed and discussed 
in relation to sonic sands [3TH33] , we place ourselves in 
the case of rapid flows, corresponding to the observation 
of a minimum velocity threshold for sound emission. Ap- 
plying discrete numerical simulation in 2D [34] [35], we 
simulate model granular flows on erodible beds; we an- 




FIG. 1: A 2D periodic flow on erodible bed tilted of an angle 
simulated by Contact Dynamics. The height of the flowing 
layer is denoted h and the length of the sumulation cell is L. 
The yellow-red shade of the grains stands for their velocity 
(color on-line); the darker-grey shade shows images of the 
flowing grains through the periodic boundary condition. 



alyze the dynamics of the simulated flows while varying 
the grains frictional properties over a range difficult to 
attain in laboratory experiments. Doing so, we show the 
existence of a spontaneous oscillation developing rapidly 
and exhibiting a well-defined frequency which increases 
like the inverse of the velocity's square root. The ampli- 
tude of the frequency is higher for lower volume fraction, 
and can thus be related to both flow velocity and grains 
friction properties. Finally, we investigate the influence 
of the spatial periodicity of the simulations on the above 
results and observe no significant effect. These results 
are discussed in relation to sonic sands. 



II. THE CONTACT DYNAMICS SIMULATIONS 

The simulations reported in this contribution were per- 
formed using the Contact Dynamics algorithm [3H [36] . 
The grains are assumed to be perfectly rigid, which trans- 
lates in a strict non-overlap condition at contact. They 
interact through a Coulombic friction law, relating the 
tangential force at each contact f t to the normal force at 
the same contact f t through the following inequality: 

\ft\ < Vfm 

where fi is the coefficient of friction at contact. In the ad- 
vent of slip motion, the equality is satisfied : \ft\ = i^fn- 
The value of \i thus controls the amount of energy dissi- 
pated by the flow through frictional contacts. Note that 
a single coefficient of friction at contact is introduced: we 
do not distinguish static and dynamical friction. In addi- 
tion, a coefficient of restitution e sets the amount of en- 
ergy dissipated in the advent of a collision, and thus con- 
trols the amount of energy dissipated by the flow through 
collisional contacts. The precise contribution of the two 
modes of dissipation (collisional and frictional) within a 
given granular flow is not straightforward to estimate, 
the multi-contact dynamics particular to dense granular 
packings being characterized by disorder and complexity 
[37] . In the following however, we will not be interested 
in the particular role of the coefficient of restitution e, 
and we will fix the value of the later to e = 0.5 coin- 
ciding with dense systems. On the contrary, our interest 
focuses on the role of the coefficient of friction \i on the 
flow dynamics. Hence, its value was alternatively set to 
fi = 0.05 (very small), fi = 0.5 (a typical value for glass 
beads is \i = 0.2 [38] while 0.5 < \i < 0.9 for singing sand 
grains [28]) and \i = 2 (high). Note that by precluding 
the use of spring-dashpot models for the contact law, the 
Contact Dynamics algorithm prevents the introduction 
of mechanical oscillators in the treatment of grains inter- 
actions, and thereby prevent artefact oscillations which 
may occur in the soft-sphere limit [33] . Details on the 
numerical method can be found in [34U36] . 
The flow configuration investigated is a 2D periodic flow 
on erodible bed tilted of an angle 6 (see Figure [I]). The 
grains show a slight size dispersity to avoid ordering ef- 
fects, too small however to induce segregation: we chose 
(dmax — dmi n ) /d = 0-4, where d is the mean grain diam- 
eter (d = 500/im). The packing is obtained by random 
rain under gravity. The flow has a periodicity in the lon- 
gitudinal (x) direction; the size L of the simulation cell 
was set to 45d in the greater part of the simulations re- 
ported here, corresponding to 3967 grains. However, a 
specific study of the effect of the value of L on the flow 
dynamics, and more specifically on the flow oscillatory 



motion, was performed and is reported in section VI 



The erodible bed condition is achieved by trapping grains 
between vertical walls at the boundary of the simulation 
cell, thus allowing the upper layer only to flow in response 
to gravity. The height of the vertical walls is fixed and 
is 30d; the height of the unconstrained (ie free to flow) 
layer is h\ depending on the volume fraction (ie depend- 
ing on the velocity), h varies between 40d and 49<i. The 




FIG. 2: a: Mean velocity of the grains u (normalized by \fgd) 
as a function of time t (normalized by \fdfg) for different 
values of the coefficient of friction at contact \i — 0.05, \i — 0.5 
and [i — 2.0, and for different tilt angles 0; b: Corresponding 
velocity profiles time-averaged over the stationary regime; the 
dotted line shows a Bagnold profile (see equation |l])); Inset: 
Velocity profiles in semi-log scale; the dotted line shows an 
exponential decay with a typical length A = 2.5gL (Color 
on-line) 



vertical position of the grains z is counted positively fol- 
lowing the upward position; the origin is set where the 
vertical walls stop (ie at the bottom of the unconstrained 
layer). For each value of /i, the tilt angle was varied so 
as to achieved different flow velocities: ranges from 16° 
to 22° for ii = 0.05, from 22° to 28° for /i = 0.5 and 
from 26° to 30° for ji = 2. In every cases, stationary 
regime is reached. The mean characteristics of the flow 
thus simulated are detailed in the next section. 



III. MEAN VELOCITY AND VELOCITY 
PROFILE 

Figure [2]- a shows the time evolution of the mean grain 
longitudinal (following x) velocity u(t) (computed over 
the total number of grains in the simulation) for the three 
values of the coefficient of friction at contact \i (namely 
0.005, 0.5 and 2) and for different values of the tilt angle 
0. In each case, stationary regime is reached after a tran- 
sient regime; we do not study the latter but focus on the 



3 



U/y/gd 



30 



20- 



10 



• u = 0.05 
a u = 0.5 
■ (x = 2.0 



~i 7T 



■ 



_ L_ 



14 16 18 20 22 24 26 28 30 32 
(deg) 

FIG. 3: Mean flow velocity U of the grains in the flowing 
layer averaged over the duration of the stationary regime, as 
a function of for \i — 0.05, \i — 0.5, and — The dotted 
lines shows the experimental relation |2]) from l a . (Color on- 
line) 



stationary regime in the following. For the corresponding 
runs, Figure [2}b shows the velocity profiles time- aver aged 
over the duration of the stationary regime. We observe 
well-developed surface flows with a velocity vanishing at 
the depth corresponding to the upper end of the vertical 
walls at the boundaries of the simulation cell. We observe 
that the velocity profiles can be reasonably approximated 
by a Bagnold scaling [39j[40](see Figure |2]-b) 



„( Z ) oc ^ (V/2 - (h - Zf' 2 ) , 



(1) 



where h is the thickness of the flowing layer, and z is 
the vertical position of the grains (counted positively in 
the upward direction). The Bagnold-like shape of u(z) 
is in contradiction with experimental and numerical ob- 
servation of linear profiles for granular flows on erodible 
beds [4TH45] . Yet the semi- log plot of the velocity profiles 
shows the existence of creep motion with an exponential 
decay over a typical length A = 2.5<i, in agreement with 
experimental findings [4TJ[42] (Figure[2]-b, inset): we con- 
clude that the erodible bed condition implemented in the 
numerical simulations reproduces the dissipative proper- 
ties of a real erodible bed configuration. However, an 
important difference exists: in real erodible bed condi- 
tions, the height of the flowing layer is selected by the 
flow itself, while in our numerical simulations, the height 
of the flowing layer is set by the height of the side walls. 
We can show that in our numerical set-up, the shape of 
the velocity profile is strongly dependent on the relative 
heights of the flowing layer and the erodible bed. This 
aspect, although of interest in connection with granular 
flow rheology, is beyond the scope of the present paper. 
We will assume that the erodible bed implemented in the 
simulations, by allowing creep motion to occur, partly re- 
produce a real erodible flow configuration. 
Averaging over the duration of the stationary regime, and 
considering only the grains in the flowing layer (ie grains 
with a positive z), we compute the mean flow velocity 
U. For the different values of the coefficient of friction 



at contact /i, U is reported as a function of the slope 
in Figure [3] The numerical flows obey the chute flow 
phenomenology observed in pQ: 



tan = tan 0\ + (tan 2 — tan 0i) exp 



f3h ygh 
~Id~U~ 



(2) 



where tan identifies with the effective frictional proper- 
ties of the flow, 0i and 2 are typical angles dependent 
on the frictional properties of the grains, i is a non di- 
mensional length scale and (3 — 0.136 (from pQ). This 
approximation is reported in Figure [3] for all values of 
the friction coefficient at contact \i. Each value of \i in- 
duces different effective frictional properties of the macro- 
scopic flow. Hence we find 1 = 14.0° and 2 = 28.7° for 
fjL = 0.05, 0i = 18.8° and 2 = 34.4° for fi = 0.5, and 
0i = 19.5° and 2 = 35.2° for /i = 2.0; we find t = 2.30, 
2.26 and 2.29 respectively. The numerical values corre- 
sponding to \i = 0.5 are quantitatively consistent with 
the experimental observation for glass beads pQ. 
Alternatively, we can show that the flow satisfy the fi(I) 
dependence by defining I = dU/(^/gh 3 ^ 2 ) [3j|5]: 



tan = tan 0\ 



tan 02 — tan 0i 



(3) 



The best fist gives J = 0.33, and 1 = 11.03° and 2 = 
33.02° for 11 = 0.05, 1 = 15.38° and 2 = 38.83° for 
fjL = 0.5, and 1 = 16.17° and 2 = 39.35° for /i = 2.0. 
The analysis of the mean velocity and velocity profiles 
of the numerical flows thus show that the latter behave 
accordingly to experimental evidence [TJ |3] . 



IV. OSCILLATORY MOTION: FREQUENCY 
AND AMPLITUDE 



The mean normal velocity V of the grains (ie in the z- 
direction), averaged over time, is expectedly zero, or the 
flow would expand infinitely and eventually be turned 
into a dilute gas. However, the instantaneous normal ve- 
locity plotted as a function of time, shows clear rapid 
oscillations around zero, revealing an oscillatory motion 
implying successions of dilation phase and compaction 
phase. Figure [4]-a shows, as an example, the time vari- 
ation of the instantaneous normal velocity of the grains 
v(t) over a short time interval for /i = 2.0 and = 30°. 
To establish whether this oscillatory dynamics involves a 
characteristic frequency, we compute the Fourier trans- 
form of v(t) over the duration of the stationary regime 
(Figure |2]-b): a well-defined peak frequency emerges; we 
denote this peak frequency / in the following. The inset 
in Figure [4[b shows the Fourier transform of the instan- 
taneous normal velocity at different depth in the flowing 
layer. We observe that the amplitude of the peak fre- 
quency rapidly decreases and eventually vanishes as we 
go deeper in the flow, whereby we conclude that the os- 
cillatory motion is a surface phenomenon rather than a 
bulk one. 

For all simulations, we measure the peak frequency / 
characterizing the oscillation over the stationary regime 
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FIG. 4: a: Instantaneous normal velocity v (normalized by 
\fgd) as a function of time t (normalized by \fdfg)\ b: Fourier 
transform of v(t) over the duration of the stationary regime 
showing a peak frequency /; Inset: Fourier transform of the 
instantaneous normal velocity at different depth in the flow. 
(Color on-line) 
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FIG. 6: a: Amplitude A (xlO 4 ) of the frequency peak as a 
function of the flow mean longitudinal velocity U (normalized 
\fgd). b: amplitude A (xlO 4 ) of the frequency peak as a 
function of the flow volume fraction <fi. (Color on-line) 
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FIG. 5: Peak frequency / (Hz) of the Fourier transform of 
the instantaneous normal flow velocity v(t) as a function of 

(^lF x 2 ( Hz )- Inset: Peak frequency / (Hz) of the 
Fourier transform of the instantaneous normal flow velocity 
v(t) as a function of the shear rate U/h (Hz). (Color on-line) 



(error bars are evaluated systematically based on the 
width of the frequency peak). We observe values ranging 
from 10 to 40 Hz, that is significantly smaller than the 
typical frequency attached to the grain size \J~gJd — 140 
Hz. The values observed for / are also below the values 
measured for sonic sands, typically of 70 to 110 Hz [29] . 
A reasonable correlation suggests that the value of / in- 
creases like i^f 1 x | ^ , although the range of values 

explored does not allow for a discussion on the shape of 
the dependence (Figure [5]): hence we will not address the 
possible origin of this correlation in the following. Inter- 
estingly, the frequency characterizing the oscillation does 
not identify with the shear rate: Figure [5] shows an in- 
verse correlation between / and U/h. 
The trend relating frequency and velocity displayed in 
Figure [5] is not dependent on the value of the coefficient 
of friction /i: the latter seems to be playing no particu- 
lar role in the value of /. However, oscillatory motion is 
not characterized only by the peak frequency /, but also 
by the amplitude A of the peak. Figure |6]-a shows the 
amplitude A (xlO 4 ) as a function of the flow mean lon- 
gitudinal velocity U: we observe a positive correlation, 
with larger values of \i inducing larger amplitudes. High 
velocities and high friction at contact thus induce larger 
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FIG. 7: Inset: Flow mean volume fraction as a function of 
the grains mean velocity U (normalized \/~gd) for the different 
values of the coefficient of friction at contact \i. Main graph: 
4> — 4> c as a function of U (normalized ^/gd), where </> c is a 
constant whose value depends on \i (see equation Q); the 
dotted line shows a linear fit with slope 0.00315. (Color on- 
line) 



oscillations. This double influence can be summed-up 
when plotting A as a function of the flow mean volume 
fraction <j> = volume-of-grains/total- volume (computed 
in the flowing layer over the duration of the stationary 
regime) (Figure [6]-b): the data collapse and show that 
larger amplitudes coincide with dilute flows, while denser 
flows induce a smaller amplitude. 



V. ON THE ROLE OF FRICTION 

The volume fraction of a granular flow is dependent 
on the dynamics: rapid flows coincide with dilute states, 
ie smaller volume fractions [3j |5j [21] . The relation be- 
tween the two involves the frictional properties of the 
material: as shown in Figure [7] (inset), plotting as a 
function of U reveals three distinct series of points co- 
inciding with each value of \i. This behavior can be de- 
scribed by the following dependence: 



<j> = c (/i) - k 



U 



gd 



(4) 



where k — 0.00315 and </> c is a decreasing function of 
fi: C = 0.86 for fi = 0.05, (j) c = 0.826 for fi = 0.5 and 
4> c = 0.800 for ji = 2 (see Figure [7]). In other words, for a 
given velocity U, higher contact friction implies smaller 
volume fraction. The dependence of <j) on the dynam- 
ics was already reported in [3j [5] in terms of the inertial 
number / = dU / (^h 3 / 2 ): (j) = <pM + (<t>m — </W)^, where 
(\>m and m are extremal values of the volume fraction. 
Relation Q thus gives a possible explanation for the role 
of contact friction in the oscillatory motion of granular 
flows: higher friction at contact between the grains is 
responsible for smaller volume fraction, thus leading to 
larger amplitude of grains vertical motion. Translated in 



a natural case, larger friction (as induced by the mineral 
coating observed at the surface of sonic sands) may lead 
to a larger amplitude of the oscillatory motion, and pos- 
sibly due to that, to the emission of an audible sound. 
In the same way, relation Q renders the fact that larger 
velocities induce more dilute flows, thus larger amplitude 
of grains vertical motion, which may explain why a min- 
imum flow velocity is necessary for sonic sands to emit 
sound [21J [23 [271 [30]. 



VI. INFLUENCE OF THE SIZE OF THE 
SIMULATION CELL 

The existence of a spatial periodicity in the numeri- 
cal systems studied here may be suspected to affect the 
time periodicity characterizing the flow dynamics. The 
length of the simulation cell being L, and the mean flow 
longitudinal velocity being J7, the obvious time scale re- 
lated to the geometrical periodicity of the system is L/U. 
If the geometrical periodicity was to affect the time pe- 
riodicity, then we would expect the latter to exhibit a 
frequency scaling like U/L. This is very different from 
the frequency / emerging from the analysis of the oscil- 
latory motion of the numerical flows, and the associated 
dependence on the velocity U (Figure [5| : / increases 

with (^ff^ x d) ' Hence, it seems unlikely that the 
periodicity of the simulation cell has a significant influ- 
ence on the oscillatory motion observed. This however 
needs clarification. 

Therefore, we perform series of simulations where L is 
varied: L = 12d, 22d, 32d, 52d, 72d, 102d, and L = 152d. 
The coefficient of friction is fixed: \i = 0.5, as well as the 
slope = 22°. As previously, we analyze the mean veloc- 
ity of the flow and the oscillatory motion visible in the 
time fluctuations of the normal velocity v(t). The peak 

frequency / is reported as a function of x | ^ 

in Figure [8]-a: although scattering occurs, the value of L 
is not significantly affecting the dependence already ob- 
served in Figure [5] 

Yet, the size of the simulation cell L might play a role 
in the amplitude of the oscillatory motion, and possi- 
bly influence its very existence. We denote A mean the 
amplitude of the Fourier transform of v(t) at high fre- 
quencies (ie 100Hz < /). The value of A mean is reported 
is Figure [8]-b as a function of L: we observe indeed that 
smaller systems favor larger oscillations. However, the 
amplitude saturates towards a minimum value for larger 
systems, showing that the oscillation is not the result of 
the system size. Accordingly, we expect the amplitude 
A of the peak frequency / to decrease and saturate for 
larger value of L. This is indeed what is observed (Figure 



There is however an aspect in the present analysis which 
prevents us from interpreting further the amplitude A of 
the oscillation: the spatial localization of the latter. In- 
deed, so far we have processed the normal velocity v(t) 
averaged over the whole flowing layer. If the oscillation 
is localized in space, part of the information is "watered 
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FIG. 8: a: Peak frequency / of the Fourier transform of 
the instantaneous normal flow velocity v(t) as a function of 

{f¥ x d) 1 2 - b: Amplitude A mean (xlO 4 ) of the Fourier 
transform of the normal velocity v(t) at high frequencies and 
the amplitude A (xlO 4 ) of the peak frequency /. 



down" by including the grains dynamics of the whole 
flow. This is all the more likely to happen that the sys- 
tem is large. It is thus probable that Figure [8]-b shows 
not only the effect of the periodicity, but also the infor- 



mation loss due to averaging over the system size. We 
can conclude nevertheless that the system size L plays 
no crucial role in the results discussed above. 



VII. CONCLUSION 

This contribution reports on numerical simulations of 
2D granular flows on erodible beds. The broad aim of 
this work is to investigate i) whether simple flows of 
model granular matter exhibits spontaneous oscillatory 
motion in generic flow conditions, and in this case, ii) 
whether the frictional properties of the contacts between 
grains may affect the existence or the characteristics of 
this oscillatory motion. The analysis of different series 
of simulations show that the flow develops an oscillatory 
motion with a well-defined frequency which increases like 
the inverse of the velocity's square root. We show that 
the oscillation is essentially a surface phenomena. The 
amplitude of the oscillation is higher for lower volume 
fractions. It can thus be related to the flow velocity, 
higher velocities favoring lower volume fraction. For the 
same reason, it is also dependent on grains friction prop- 
erties: indeed, large contact friction is found to induce 
lower volume fraction, and thus larger amplitude. The 
study of the influence of the periodic geometry of the 
simulation cell shows no significant effect. 
Although one should be careful while drawing analogy 
between simple numerical models and nature, it is inter- 
esting to discuss these results in relation to sonic sands. 
Indeed, this work suggests that surface oscillation is likely 
to develop during the flow of granular matter, and that 
its amplitude is dependent on both velocity and grains 
properties, in agreement with observation. How this os- 
cillation develops into a loud audible sound is beyond 
the reach of simple 2D discrete simulations; however they 
provide an interesting insight into basic mechanisms that 
may be relevant to the complex question of how sand may 
produce sound. 
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